home *** CD-ROM | disk | FTP | other *** search
/ Games of Daze / Infomagic - Games of Daze (Summer 1995) (Disc 1 of 2).iso / x2ftp / msdos / math / fft142 / fft.c < prev    next >
C/C++ Source or Header  |  1991-06-06  |  19KB  |  697 lines
  1. /*    name...
  2.         fft
  3.  
  4.     bugs...
  5.         Should do transforms on several functions (-b option).
  6.         use large model.
  7.         eliminate use of extra few entries, permitting use of exactly
  8.             64K arrays (4 bytes * 8K entries * 2 arrays)
  9.  
  10.     performance...
  11.         4096 point transform in 1:13.71 with 6 MHz 80286 and 4 MHz 80287.
  12.  
  13.     bugs...
  14.         check for unequally spaced data is too strict.
  15.  
  16.  
  17.     history...
  18.         6 Jun 91    Ver 1.42: can read data from stdin.  Uses args.c
  19.                     to read switches and echo command line.
  20.         19 Apr 89    Ver 1.41: forcing space between columns in output.
  21.         1 Jun 87    Checking for unequally spaced data.
  22.         -- 1.40 --
  23.         11 May 87    -o option installed, -p    option accepts total number
  24.                     as well as padding factor.
  25.         8 May 87    Reading into doubles and performing range check
  26.                     before converting to floats.
  27.         1 May 87    calculating with floats rather than doubles.
  28.                     performing 4096 point transforms.
  29.                     -a option suppresses frequencies on output
  30.         -- 1.3 --
  31.         29 Apr 87    Fixed padding factor calculation
  32.                     added -m option to print only freq and mag
  33.                     added -z option to shift origin
  34.                     enhanced help display
  35.         6 Apr 87    Default edge weight .1 instead of .9,
  36.                     Default abscissa step 1 (instead of 0!).
  37.         20 Mar 87    Conversion to/from dBs or to magnitude/phase removed.
  38.                     Default output to STDOUT.
  39.                     Removed prompts.
  40.         8 Feb 87    Inverse FFT code split off into separate file.
  41.         7 Nov 86    -o omits phase printout, -p specifies padding factor.
  42.             padding factor defaults to 1/4 of maximum, # points defaults
  43.             to all.
  44.         6 Nov 86    Switch s sets up windowing.
  45.         2 Nov 86    Ported to DeSmet C, input & output files in ASCII,
  46.             options read from command line.
  47.         12 Jul 84  Printing time between samples.
  48.         11 Jul 84  One sided windows, degree 6, 10, or 16.
  49.             Calculating energy before & after windowing.
  50.         10 Jul 84  No longer printing out y values. Announcing
  51.             program version.  One sided windowing, degree 16.
  52.             Allows peak to be any value from 200 to 22000 without
  53.             rescaling.  No longer discarding zero frequency point.
  54.         26 Jun 84  Weight at window edge (epsilon) specified by
  55.             user.  Number of input data points to be retained
  56.             is specified by user.
  57.         25 Jun 84  Supergaussian windowing (degree 6) is optional,
  58.             beeping after finishing transform.
  59.         22 Jun 84  Scaling real & imaginary parts to the range
  60.             10900...22000.  Padding by a factor of 4, then 
  61.             discarding 50% of the data (high frequency values) so
  62.             512 magitude values are retained.  No "magnitude"
  63.             values printed out (not even zeros).
  64.         18 Jun 84  writing out first magnitude, then phase.
  65.             Correctly calculating time step.
  66.         15 Jun 84  setting phase to zero.
  67.         12 Jun 84  scaling so largest element is 2048.
  68.             Calculating magnitude & phase.
  69.         May 84     written
  70. */
  71. #include <stdio.h>
  72. #include <math.h>
  73.  
  74. #define VERSION "1.42"
  75.  
  76.  
  77.  
  78. #define BIGFLOAT 6.e38    /* no floats larger than this  */
  79. #define ENTRIES 4098
  80. #define MAXLABELS 50
  81. #define BUFSIZE 200
  82. #define DEFAULT_EDGE_WEIGHT .1
  83. #define FLOAT float        /* change to "double" to restore higher precision */
  84.  
  85. FILE *ifile=stdin, *ofile=stdout;
  86.  
  87. char buffer[128], iname[35], oname[35];
  88. char buf[BUFSIZE];
  89. char *label_array[MAXLABELS];
  90. char **p_text=label_array;
  91.  
  92.  
  93. int m,            /* number of input values to be retained */
  94. n,                /* number of data points to be Fourier transformed */
  95. *nnn,            /* points to n */
  96. breaking=0,        /* nonzero if finding separate transforms for several functions */
  97. labels,            /* number of labels stored */
  98. super=0,        /* if nonzero, the degree of supergaussian window desired */
  99. automatic_abscissas, /* nonzero if abscissas to be calculated */
  100. abscissa_arguments=0,
  101. shifting=0,    /* nonzero if data to be shifted */
  102. pad_factor=0,    /* nonzero if specific padding factor was requested */
  103. keeping=0,        /* if nonzero, the number of input values to be kept */
  104. magnitudes=0,    /* nonzero if only magnitudes are to be written */
  105. f_arguments=0,
  106. x_arguments=0;
  107.  
  108. int standard_input=0;        /* nonzero if data is coming from standard input */
  109. int index_array[MAXLABELS];    /* indexes into x and y */
  110. int *p_data=index_array;
  111.  
  112. FLOAT *x,        /* time or frequency values */
  113. *y;                /* voltages (on input) or magnitudes (on output) */
  114. double abscissa,
  115. abscissa_step=1.,
  116. before, after,    /* energy in data before & after windowing */
  117. origin,            /* abscissa value to be treated as time zero */
  118. wt=0.,            /* weight on last data point kept */
  119. output=0.,        /* if nonzero, the number of calculated values to be printed */
  120. fmin, fmax,
  121. xmin, xmax;        /* first & last data points to be used */
  122.  
  123. double energy(), atan2();
  124.  
  125. main(argc, argv) int argc; char **argv;
  126. {    int i, nn, windowing;
  127.  
  128.     double amp, phase, factor, q, pi, time, dt;
  129.     double this, big, sweep, scale, yr, yi;
  130.  
  131.     {double junk;
  132.     sscanf("3.4","%lf",&junk);        /* workaround for DESMET sscanf bug */
  133.     }
  134.  
  135.     argc = args(argc, argv);        /* fetch switches */
  136.  
  137.     read_data(argc, argv);            /* read input data */
  138.  
  139. /*    puts(" time function\n"); listt(y, 1, 10); */
  140.     nn=n;
  141.         {double dt;
  142.         dt=x[2]-x[1];
  143.         if(nn*fabs(x[nn]-x[nn-1]-dt) + fabs(x[nn]-x[1]-(nn-1)*dt) > .001*nn*dt)
  144.             {printf("fft: times not equally spaced\n");
  145.             printf("x[1]=%f, x[2]=%f, x[%d]=%f, x[%d]=%f, dt=%f",
  146.                 x[1], x[2], nn-1, x[nn-1], nn, x[nn], dt);
  147.             exit();
  148.             }
  149.         }
  150.     nn=pad();
  151. /*    puts(" time function after padding\n"); listt(y, 1, 10); */
  152.     before=energy();
  153.     windowing = wind16() || wind10() || wind6();
  154.     if(windowing)
  155.         {after=energy();
  156.         fprintf(stderr, "Energy loss due to windowing = %5.2f%% \n",
  157.             100.*(before-after)/before);
  158.         }
  159.     else after=before;
  160.     fprintf(stderr, "Energy (sum of squares of data%s) %s= %10.3g \n", 
  161.     (!automatic_abscissas || abscissa_arguments) ? ", times time step" : "",
  162.     windowing?"after windowing ":"", after);
  163. /*    puts(" time function\n"); listt(y, 1, 10); listt(y, nn-10, nn); */
  164.     if(shifting) shift(nn);
  165.     
  166.     fprintf(stderr, "computing FFT...\n");
  167.     fft(&nn, x, y);
  168.  
  169.     if(output)
  170.         {if(output<32) output = nn/output;
  171.         if(output<nn) nn=output;
  172.         }
  173.     if(abscissa_arguments)
  174.         fprintf(ofile, "; output frequency step is %15.8g", x[2]-x[1]);
  175.     i=1;
  176.     while(1)
  177.         {if(!automatic_abscissas)  fprintf(ofile, "%15.6g ", x[i]);
  178.         yr=y[2*i-1]; yi=y[2*i];
  179.         if(magnitudes)  fprintf(ofile, "%15.6g", sqrt(yr*yr + yi*yi));
  180.         else            fprintf(ofile, "%15.6g %15.6g", yr, yi);
  181.         if(++i>nn) break;
  182.         fprintf(ofile, "\n");
  183.         }
  184.     fprintf(ofile, " \"\"\n");
  185.     if(ofile!=stdout)
  186.         {fclose(ofile);
  187.         }
  188. }
  189.  
  190. double energy()
  191. {    double v, e;
  192.     int i;
  193.     i=1;
  194.     e=0.;
  195.     while(i<=m)
  196.         {v=y[i++]; e=e+v*v;
  197.         }
  198.     return e*(x[2]-x[1]);
  199. }
  200.  
  201. pad()        /* pad the data with zeros */
  202. {    int i, k, num, max;
  203.     char c;
  204.     if(keeping)
  205.         m=keeping;
  206.     else 
  207.         m=n;
  208. /*
  209.         while(1)
  210.             {fprintf(stderr,
  211.             "Last point to include (2...%d, default %d) ? ", n, n);
  212.             gets(buf);
  213.             if(buf[0]) sscanf(buf, "%d", &m);
  214.             else m=n;
  215.             if((m>=2) && (m<=n)) break;
  216.             }
  217. */
  218.     if(pad_factor)
  219.         {if(pad_factor>=32) num=(pad_factor/10)*7;
  220.         else if(pad_factor>0) num=((m*pad_factor)/10)*7;
  221.         else num=m*4;    /*  default is a padding factor between 4 and 8 */
  222.         }
  223.     else num=m*4;    /*  default is a padding factor between 4 and 8 */
  224.     for (k=1; k<num; k <<= 1) {}    /* find next power of 2 */
  225.     while (k>ENTRIES) {k >>= 1;}    /* back off if too many for buffer */
  226.     fprintf(stderr, "transforming %d points, for actual padding factor of %3.1f\n", 
  227.         k, k/(double)m);
  228.     i=m;
  229.     while(i<k)    /* will have k values after this */
  230.         {y[++i]=0.;}
  231.     return k;
  232. }
  233.  
  234. shift(n) int n;        /* shift data so time "origin" is at the beginning of the array */
  235. {    int k;
  236.     k=(int)((origin-x[1])/(x[2]-x[1]) + 1.5);
  237.     if(k<1 || k>n) 
  238.         {fprintf(stderr, "specified origin isn't within abscissa range");
  239.         exit(1);
  240.         }
  241.             /* need to swap y[1]...y[k-1] with y[k]...y[n] */
  242.     reverse(1, k-1);
  243.     reverse(k, n);
  244.     reverse(1, n);
  245. }
  246.  
  247. reverse(i, j) int i,j;  /* reverse y[i]...y[j]  */
  248. {    double t;
  249.     while(i<j) {t=y[i]; y[i]=y[j]; y[j]=t; i++; j--;}
  250. }
  251.  
  252. wind16()
  253. {    char c;
  254.     int j;
  255.     double p, p2, dp, scale, epsilon;
  256.  
  257. /*    if(!super)
  258.         {while(1)
  259.             {puts("Perform one-sided supergaussian windowing, degree 16? ");
  260.             c=toupper(getchar()); putchar('\n');
  261.             if(c=='N')return 0;
  262.             if(c=='Y')break;
  263.             }
  264.         }
  265.     else
  266. */
  267.          if(super!=16) return 0;        
  268. /*
  269.     if(wt==0.)
  270.         {fprintf(stderr, "Weight at window edge (0 < epsilon < 1, default=%f) ? ",
  271.             DEFAULT_EDGE_WEIGHT);
  272.         gets(buffer);
  273.         if(buffer[0]) epsilon=atof(buffer); else epsilon=DEFAULT_EDGE_WEIGHT;
  274.         }
  275.     else 
  276. */        epsilon=wt;
  277.     p=exp(log(-log(epsilon))/16.);    /* (-ln(epsilon))**(1/16) */
  278. /*    fprintf(stderr, "ratio C/A = %7.4f \n", p);  */
  279.     dp=p/(m-1);
  280.     j=m;
  281.     while(j)
  282.         {p2=p*p; p2=p2*p2; p2=p2*p2;    /* p**8 */
  283.         y[j] *= exp(-p2*p2);    /* p**16 */
  284.         p -= dp; --j;
  285.         }
  286.     return 1;
  287. }
  288.  
  289. wind10()
  290. {    char c;
  291.     int j;
  292.     double p, p2, p4, dp, scale, epsilon;
  293.  
  294. /*    if(!super)
  295.         {while(1)
  296.             {puts("Perform one-sided supergaussian windowing, degree 10? ");
  297.             c=toupper(getchar()); putchar('\n');
  298.             if(c=='N')return 0;
  299.             if(c=='Y')break;
  300.             }
  301.         }
  302.     else
  303. */
  304.          if(super!=10) return 0;        
  305. /*
  306.     if(wt==0.)
  307.         {fprintf(stderr,"Weight at window edge (0 < epsilon < 1, default=%f) ? ",
  308.             DEFAULT_EDGE_WEIGHT);
  309.         gets(buffer);
  310.         if(buffer[0]) epsilon=atof(buffer); else epsilon=DEFAULT_EDGE_WEIGHT;
  311.         }
  312.     else 
  313. */
  314.         epsilon=wt;
  315.     p=exp(log(-log(epsilon))/10.);    /* (-ln(epsilon))**(1/10) */
  316. /*    fprintf(stderr, "ratio C/A = %7.4f \n", p);  */
  317.     dp=p/(m-1);
  318.     j=m;
  319.     while(j)
  320.         {p2=p*p; p4=p2*p2;    /* p**2, p**4 */
  321.         y[j] *= exp(-p2*p4*p4);    /* p**10 */
  322.         p -= dp; --j;
  323.         }
  324.     return 1;
  325. }
  326.  
  327. wind6()
  328. {    char c;
  329.     int j;
  330.     double p, p2, dp, scale, epsilon;
  331.  
  332. /*    if(!super)
  333.         {while(1)
  334.             {puts("Perform one-sided supergaussian windowing, degree 6? ");
  335.             c=toupper(getchar()); putchar('\n');
  336.             if(c=='N')return 0;
  337.             if(c=='Y')break;
  338.             }
  339.         }
  340.     else
  341. */
  342.          if(super!=6) return 0;        
  343. /*
  344.     if(wt==0.)
  345.         {fprintf(stderr,"Weight at window edge (0 < epsilon < 1, default=%f) ? ",
  346.             DEFAULT_EDGE_WEIGHT);
  347.         gets(buffer);
  348.         if(buffer[0]) epsilon=atof(buffer); else epsilon=DEFAULT_EDGE_WEIGHT;
  349.         }
  350.     else 
  351. */
  352.         epsilon=wt;
  353.     p=exp(log(-log(epsilon))/6.);    /* (-ln(epsilon))**(1/6) */
  354. /*    fprintf(stderr, "ratio C/A = %7.4f \n", p);  */
  355.     dp=p/(m-1);
  356.     j=m;
  357.     while(j)
  358.         {p2=p*p*p;  /* p**3 */
  359.         y[j] *= exp(-p2*p2);    /* p**6 */
  360.         p -= dp; --j;
  361.         }
  362.     return 1;
  363. }
  364.  
  365. listt(y, i, j) FLOAT y[]; int i, j;    /* display y[i] through y[j]  */
  366. {    while (i<=j) {printf("%4d %15.9f \n", i, y[i]); ++i;}
  367. }
  368.  
  369.  
  370. read_data(argc, argv) int argc; char **argv;
  371. {    int i, j, length;
  372.     double xx, yy, d, *pd, sum;
  373.     char *s, *t, *stop;
  374.     char *malloc();
  375.     char *strchr();
  376.  
  377.     x=(FLOAT *)malloc(ENTRIES*sizeof(FLOAT));
  378.     y=(FLOAT *)malloc(ENTRIES*sizeof(FLOAT));
  379.     if(x==0 || y==0) {fprintf(stderr, "can\'t allocate buffer"); exit(1);}
  380.     argc--; argv++;
  381.     if(strchr(argv[0], '?')) help();
  382. /*            open input file         */
  383.     if(argc && *argv[0]!='-')
  384.         {ifile=fopen(argv[0], "r");
  385.         if(ifile==0) {fprintf(stderr, "file %s not found\n", argv[0]); exit(1);}
  386.         strcpy(iname, argv[0]);
  387.         argc--; argv++;
  388.         }
  389. /*            open output file         */
  390.     if(argc && *argv[0]!='-')
  391.         {strcpy(oname, argv[0]);
  392.         argc--; argv++;
  393.         unlink(oname);
  394.         if((ofile=fopen(oname, "w"))==0)
  395.             {fprintf(stderr, "can\'t open output file %s", oname);
  396.             exit(1);
  397.             }        
  398.         }
  399.     p_data[0]=-1;
  400.  
  401.     fprint_cmd(ofile, "; fft %s\n");
  402.  
  403.     i=1;        /* note x[0] and y[0] aren't defined */
  404.     while(i<ENTRIES)
  405.         {if(fgets(buf, BUFSIZE, ifile)==0)
  406.             {if(standard_input) {fclose(ifile); ifile=fopen("/dev/con", "r");}
  407.             break;
  408.             }
  409.         t=buf;
  410.         while(*t && isspace(*t)) t++;
  411.         if(*t == 0) continue;        /* skip blank lines */
  412.         buf[strlen(buf)-1]=0;        /* zero out the line feed */
  413.         if(buf[0]==';')                /* skip comment */
  414.             {fprintf(ofile, "%s\n", buf); continue;
  415.             }
  416.         if(t = strchr(buf, ';')) *t=0;    /* zap same-line comment */
  417.         if(automatic_abscissas)
  418.             {xx=abscissa;
  419.             abscissa+=abscissa_step;
  420.             sscanf(buf, "%lf", &yy); 
  421.             }
  422.         else
  423.             {
  424.             sscanf(buf, "%lf %lf", &xx, &yy);
  425.             }
  426.         range_check(xx); range_check(yy);
  427.         x[i]=xx; y[i]=yy;        /* convert doubles to floats here */
  428.         s=buf;                      /* start looking for label */
  429.         while(*s==' ')s++;            /* skip first number */
  430.         while(*s && (*s!=' '))s++;
  431.         if(!automatic_abscissas)    /* skip second number */
  432.             {while(*s==' ')s++;
  433.             while(*s && (*s!=' '))s++;
  434.             }
  435.         while(*s==' ')s++;
  436.         if((length=strlen(s))&&(labels<MAXLABELS))
  437.             {if(*s=='\"')           /* label in quotes */
  438.                 {t=s+1;
  439.                 while(*t && (*t!='\"')) t++;
  440.                 t++;
  441.                 }
  442.             else                    /* label without quotes */
  443.                 {t=s;
  444.                 while(*t && (*t!=' '))t++;
  445.                 }
  446.             *t=0;
  447.             length=t-s;
  448.             p_data[labels]=i;
  449.             p_text[labels]=(char *)malloc(length+1);
  450.             if(p_text[labels]) strcpy(p_text[labels++], s);
  451.             }
  452.         i++;
  453.         }
  454.     n=i-1;
  455. /*    if(ifile==stdin) {close(ifile); ifile=fopen();} */
  456.     if(breaking && (!labels || p_data[labels-1]!=n-1))
  457.         {p_data[labels]=i-1;
  458.         if(p_text[labels]=(char *)malloc(1)) *p_text[labels++]=0;
  459.         }
  460. }
  461.  
  462. /* check whether number is too big for a float */
  463. range_check(x) double x;
  464. {    if(fabs(x)>BIGFLOAT)
  465.         {fprintf(stderr, "input number too large: %f\n", x);
  466.         exit(1);
  467.         }
  468. }
  469.  
  470.  
  471. int streq(a,b) char *a,*b;
  472. {    while(*a)
  473.         {if(*a!=*b)return 0;
  474.         a++; b++;
  475.         }
  476.     return (*b==0);
  477. }
  478.  
  479. /* get_parameter - process one command line option
  480.         (return # parameters used) */
  481. get_parameter(argc, argv) int argc; char **argv;
  482. {    int i;
  483.     if(streq(*argv, "-a"))
  484.         {i=get_double(argc, argv, 1, &abscissa_step, &abscissa, &abscissa);
  485.         abscissa_arguments=i-1;
  486.         automatic_abscissas=1;
  487.         return i;
  488.         }
  489.     else if(streq(*argv, "-b")) {breaking=1; return 1;}
  490.     else if(streq(*argv, "-s"))
  491.         {super=10; wt=DEFAULT_EDGE_WEIGHT;
  492.         if((argc>1) && numeric(argv[1])) 
  493.             {super=atoi(argv[1]);
  494.             if((super!=16) && (super!=10) && (super!=6))
  495.                 {fprintf(stderr, "windowing only of degrees 6, 10, & 16 available");
  496.                 exit(1);
  497.                 }
  498.             if((argc>2) && numeric(argv[2]))
  499.                 {wt=atof(argv[2]);
  500.                 if((wt<=0.)||(wt>=1.))
  501.                     {fprintf(stderr, "wt=%f, need  0 < wt < 1", wt); exit(1);
  502.                     }
  503.                 return 3;
  504.                 }
  505.             return 2;
  506.             }
  507.         return 1;
  508.         }
  509.     else if(streq(*argv, "-n"))
  510.         {if((argc>1) && numeric(argv[1])) keeping=atoi(argv[1]);
  511.         if(keeping<2 || keeping>ENTRIES)
  512.             {fprintf(stderr, "switch -n option out of range 2...%d\n",ENTRIES);
  513.             exit(1);
  514.             }
  515.         return 2;
  516.         }
  517.     else if(streq(*argv, "-o"))
  518.         {if((argc>1) && numeric(argv[1])) output=atof(argv[1]);
  519.         if(output<1 || output>ENTRIES)
  520.             {fprintf(stderr, "switch -o option out of range 1...%d\n",ENTRIES);
  521.             exit(1);
  522.             }
  523.         return 2;
  524.         }
  525.     else if(streq(*argv, "-p"))
  526.         {if((argc>1) && numeric(argv[1])) pad_factor=atoi(argv[1]);
  527.         return 2;
  528.         }
  529.     else if(streq(*argv, "-f"))
  530.         {i=get_double(argc, argv, 2, &fmin, &fmax, &fmax);
  531.         f_arguments=i-1;
  532.         return i;
  533.         }
  534.     else if(streq(*argv, "-t"))
  535.         {i=get_double(argc, argv, 2, &xmin, &xmax, &xmax);
  536.         x_arguments=i-1;
  537.         return i;
  538.         }
  539.     else if(streq(*argv, "-m"))
  540.         {magnitudes=1;
  541.         return 1;
  542.         }
  543.     else if(streq(*argv, "-z"))
  544.         {origin=0.;
  545.         i=get_double(argc, argv, 1, &origin, &fmax, &fmax);
  546.         shifting=1;
  547.         return i;
  548.         }
  549.     else gripe(argv);
  550. }
  551.  
  552. char *message[]={
  553. "fft   version ", VERSION,
  554. " - calculate fast Fourier transform of time domain curve\n",
  555. "usage: fft  [infile  [outfile]]  [options]\n",
  556. "options are:\n",
  557. "  -a  [step]       automatic abscissas \n",
  558. /*
  559. "  -b               break input after each label,  \n",
  560. "                   find separate transforms\n",
  561. "  -f  min [max]    minimum and maximum frequencies\n",
  562. "  -t  min [max]    minimum and maximum times\n",
  563. */
  564. "  -m               print only frequencies & magnitudes\n",
  565. "  -p  num          pad data by factor of num\n",
  566. "  -n  num          keep only first  num  input data points\n",
  567. "  -o  num          print out num values\n",
  568. "  -s  [deg [wt]]   perform supergaussian windowing of\n",
  569. "             degree deg (6, 10, or 16, default 10) with\n",
  570. "             weight wt on last point (default .9)\n",
  571. "  -z  origin       subtract  origin  from each abscissa value\n",
  572. 0
  573. };
  574.  
  575. /*    name...
  576.         fft
  577.  
  578.     purpose...
  579.         perform forward or inverse FFT
  580.  
  581.     history...
  582.         May 84  translated from FORTRAN
  583.         12 Jun 84  diagnostic printouts shortened
  584. */
  585.  
  586.     int ip, n2, i, j, k, nu, nn;
  587.     int k1n2, l, nu1;
  588.     double arg, p, q, s, c;
  589.     double ninv, time, dt, freq, df;
  590.     double xplus, xminus, yplus, yminus;
  591.     double uplus, uminus, vplus, vminus;
  592.  
  593.     /* format time domain function for calculation
  594.         of its fast Fourier transform */
  595. /*
  596. on input:
  597. n = a power of 2
  598. x[1...n] has  t(0)  t(1)  t(2)  ...  t(n-1)
  599. y[1...n] has  q(0)  q(1)  q(2)  ...  q(n-1)
  600.  
  601. on output:
  602. n = (n/2+1), one more than a power of 2
  603. x[1...n]  has  f(0)  f(1)  f(2)  ...  f(n-1)
  604. y[1...2n] has  Qr(0) Qi(0)   Qr(1) Qi(1)   Qr(2) Qi(2)  ...  Qr(n-1)  Qi(n-1)
  605. */
  606. fft(nnn, x, y) int *nnn; FLOAT x[], y[];
  607. {    n= *nnn;
  608.     if(n<=0){puts("fft: illegal # points"); return;}
  609.     dt=x[2]-x[1];
  610.     df=1./(n*dt);
  611.     n2=4; nu=0;
  612.     while(++nu<=20)
  613.         {if(n2==n)break;
  614.         n2=n2+n2;
  615.         }
  616.     if(nu>20){puts("fft: n not a power of 2"); return;}
  617. /*    puts("fft: n is ok, nu= %d\n", nu);    */
  618.     n=n/2; i=0;
  619.     while(++i<=n) {x[i]=y[2*i]; y[i]=y[2*i-1];}
  620. /*
  621.     printf("fft: data reformatted\n");
  622.     i=0; 
  623.     while(++i<=10) {printf("%6d %15.6f %15.6f\n", i, y[i], x[i]);}
  624. */
  625.     fft2(y, x, n, nu);
  626. /*
  627.     puts("fft: transform calculated\n"); 
  628.     i=0; 
  629.     while(++i<=10) {printf("%6d %15.6f %15.6f\n", i, y[i], x[i]);}
  630. */
  631.     ip=n+1; y[ip]=y[1]; x[ip]=x[1];
  632.     ninv=3.14159265358979/n;
  633.     n=n/2; n2=n+1; i=0;
  634.     while(++i<=n2)
  635.         {arg=ninv*(i-1); s=sin(arg); c=cos(arg);
  636.         uplus=y[i]+y[ip]; uminus=y[i]-y[ip];
  637.         vplus=x[i]+x[ip]; vminus=x[i]-x[ip];
  638.         p= c*vplus-s*uminus; q= -s*vplus-c*uminus;
  639.         y[i]=.5*(uplus+p); y[ip]=.5*(uplus-p);
  640.         x[i]=.5*(q+vminus); x[ip]=.5*(q-vminus);
  641.         --ip;
  642.         }
  643.     n=n+n+1; ip=n;
  644.     while(ip) {y[ip+ip-1]=dt*y[ip]; --ip;}
  645.     freq=0.; i=0;
  646.     while(++i<=n) {y[i+i]=x[i]*dt; x[i]=freq; freq=freq+df;}
  647.     *nnn=n;
  648. }
  649.  
  650. fft2(xreal, ximag, n, nu)
  651. FLOAT xreal[], ximag[];
  652. int n, nu;
  653. {    double treal, timag;
  654.     FLOAT *xr1, *xr2, *xi1, *xi2;
  655.  
  656.     n2=n/2; nu1=nu-1; k=0; l=0;
  657.     while(++l<=nu)
  658.         {while(k<n)
  659.             {p=ibitr(k>>nu1, nu);
  660.             arg=3.14159265358979*2*p/n;
  661.             c=cos(arg); s=sin(arg); i=0;
  662.             while(++i<=n2)
  663.                 {++k; k1n2=k+n2;
  664.                 /* initialize four pointers */
  665.                 xr1=xreal+k; xr2=xreal+k1n2;
  666.                 xi1=ximag+k; xi2=ximag+k1n2;
  667.                 treal= *xr2*c+ *xi2*s;
  668.                 timag= *xi2*c- *xr2*s;
  669.                 *xr2= *xr1-treal;
  670.                 *xi2= *xi1-timag;
  671.                 *xr1= *xr1+treal;
  672.                 *xi1= *xi1+timag;
  673.                 }
  674.             k=k+n2;
  675.             }
  676.         k=0; --nu1; n2=n2/2;
  677.         }
  678.     k=0;
  679.     while(++k<=n)
  680.         {i=ibitr(k-1, nu)+1;
  681.         if(i>k)
  682.             {treal=xreal[k];   timag=ximag[k];
  683.             xreal[k]=xreal[i]; ximag[k]=ximag[i];
  684.             xreal[i]=treal;    ximag[i]=timag;
  685.             }
  686.         }
  687. }
  688.  
  689. /* reverse the last nu bits of j */
  690. ibitr(j, nu) int j, nu;
  691. {    int ib;
  692.     ib=0;
  693.     while(nu--){ib=(ib<<1)+(j&1); j=j>>1;}
  694.     return ib;
  695. }
  696.  
  697.